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We numerically study the three-dimensional generalization of the kinetically constrained East 
model, the North-or-East-or-Front (NEF) model. We characterize the equilibrium behaviour of the 
NEF model in detail, measuring the temperature dependence of several quantities: a-relaxation 
time, distributions of relaxation times, dynamic susceptibility, dynamic correlation length, and 
four-point susceptibility. We show that the NEF model describes quantitatively experimental ob- 
servations over an exceptionally wide range of timescales. We illustrate this by fitting experimental 
data obtained both in the mildly supercooled regime by optical Kerr effect, and close to the glass 
transition by dielectric spectroscopy. 
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I. INTRODUCTION 

Supercooled liquids can be studied using theoretical 
models where the density field dynamics is replaced 
by a coarse-grained mobility field evolving with simple 
empirical rules based on the idea of dynamic facilita- 
tion [jj]. Kinetically constrained models (KCM), such as 
the Fredrickson-Andersen (FA) model |2| and the East 
model 0, , have been shown to reproduce well the phe- 
nomenology of supercooled liquids. An important feature 
of these models is that their dynamics becomes spatially 
correlated 0, |(|, i.e., growing timescales are accompa- 
nied by growing dynamical lengthscales, giving rise to 
dynamic heterogeneity as observed in experiments and 
simulations 

In this paper, we study by means of extensive numer- 
ical simulations the three-dimensional generalization of 
the East model 0, the North-or-East-or-Front (NEF) 
model [llj. This model is a realization of the 'ar- 
row' model of Ref. [l^, but with an externally im- 
posed directional preference while in the arrow model 
excitations carry an orientation which locally determines 
the directionality of the dynamics, making the system 
isotropic 12J. These models were shown to reproduce 
the dynamic behaviour of a wide range of fragile materi- 
als Irj. The present study is complementary to that of 
Refs. [l3l Il4| that considered systems with fully isotropic 
dynamic facilitation, such as the three-dimensional FA 
model, which describes strong materials. 

The principal aims of this work are the following: 

(i) Quantitative characterization of the fragile limit. 
We report numerical results for a wide range of dynamic 
observables in the NEF model, describing in detail distri- 
butions of the relevant timescales and lengthscales. We 
also confirm some of the predictions for the arrow model 
of Ref. 0. 

(ii) Comparison to the strong limit. We contrast 
our results for the NEF model to those obtained in 
Refs. [liillj] which dealt with the strong limit of isotropic 
facilitation, such as the three-dimensional FA model. 

(Hi) Comparison to experimental data. We use the 



NEF model results to fit experimental data over a wide 
range of relaxation times, covering both the onset of su- 
percooling, and the regime close to the glass transition. 

The paper is organized as follows. We first describe 
the model and some technical details in Sec.[n] We then 
turn to the study of the relevant timescales, how they are 
distributed and evolve with temperature in Sec. IIIII We 
study spatial aspects of the dynamic in Sec. II VI dynamic 
heterogeneity and four-point correlations. We compare 
our numerical results to experimental data in Sec.0 and 
conclude the paper in Sec. IVII 



II. MODEL AND SIMULATION DETAILS 

The NEF model is defined by the Hamiltonian 

N 

H = ^2m, (i) 

i=l 

where rij = 0,1 are N — L 3 binary variables defined 
on each site of a cubic lattice of linear size L, which 
has periodic boundary conditions. Physically, = 1 
(rij = 0) describes at a coarse-grained level a site which 
is mobile (immobile). 

The dynamics of the model can be written as follows, 

Cj c 

n t = . m = 1, (2) 

Ci (l-c) 

where 

c = (n * ) = i + ex P (i/ry (3) 

is the mean concentration of mobile regions and Ci is the 
kinetic constraint, 

Ci = l- U(l-n 3 ). (4) 

(hi)' 
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The prime indicates that the product is restricted to 
those three sites, j, which are north, front, and east near- 
est neighbours of site i. The kinetic constraint (0J implies 
that site i cannot change state unless at least one of those 
three neighbours is in the state rij = 1. The system is 
therefore more constrained than the corresponding one- 
spin facilitated FA model, where any of its six neighbours 
can facilitate the site i. In the arrow model, anisotropy is 
only local and results physically from the empirical obser- 
vation that facilitation has locally a preferred direction 
that may persist over a certain length scale even though 
the system is on average totally isotropic ^1 • The NEF 
model corresponds therefore to having an infinitely large 
persistence length of dynamic facilitation. This is physi- 
cally unrealistic, but constitutes a useful limit effectively 
describing systems with a large but finite persistence of 
dynamic facilitation [Tlj . 

We performed Monte-Carlo simulations of the NEF 
model for several temperatures in the range T £ 
[0.15, 5.0], covering over eleven orders of magnitude in re- 
laxation times. Such large time scales can be simulated 
using a continuous time algorithm 

E3. As we shall see, 

this dynamical slow-down is accompanied by the growth 
of dynamic spatial correlations, and one has to pay at- 
tention to possible finite size effects . This is however 
not a serious problem in the present model because even 
at very large time scales spatial correlations are not very 
large (see below) . At the lowest temperature simulated a 
system size of TV = 32 3 proved to be large enough. Com- 
pare this to the system size N = 160 3 used in the strong 
limit to simulate similar timescales 0] . Equilibrium be- 
haviour is trivial due to the non-interacting Hamiltonian 
(Q. It is straightforward to produce independent equili- 
brated initial configurations. Simulations therefore only 
consist of production runs. Averages are then performed 
over truly independent initial configurations. 

III. TIME SCALES 

A. Global dynamics 

We first consider the spatially averaged dynamics, 
which may be probed via the mean persistence function, 



where (i) is the single-site persistence function at time 
t, which takes the value 1 if site i has not flipped up 
to time t, and the value otherwise. As discussed in 
Ref. [ljj, the persistence function P(t) can be seen as 
the analog of the self-intermediate scattering function, 
F s (k,t), of a supercooled liquid at wavelengths compa- 
rable to the particle diameter |2(j. Fig. ^ shows, as ex- 
pected, that the dynamics slows down markedly when 
temperature is decreased below T a « 1.0, indicating the 
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FIG. 1: Persistence function, Eq. 10, for temperatures de- 
creasing from left to right: T = 4.0, 2.0, 1.0, 0.6, 0.4, 0.3, 
0.23, 0.2, 0.18, 0.16, and 0.15. 

onset of slow dynamics in this model 0, • Physically, 
T Q corresponds to the energy scale of the problem, i.e. 
the energy needed to create an excitation, see Eq. ©. 

The shape of the mean persistence function will be 
discussed below in some detail. For the moment, we 
simply note that the long-time decay is non-exponential, 
as is commonly observed in supercooled liquids. More- 
over, the short-time dynamics also presents non-trivial 
features, see for instance the curve for T = 0.15 in Fig.^ 

The mean relaxation time, t(T), is extracted via the 
usual relation P(t) — e _1 . In the present context, r rep- 
resents also the a-relaxation of the model. The tempera- 
ture dependence of r is shown in Fig. [5] The main obser- 
vation from this figure is that the mean-relaxation time 
grows with decreasing temperature in a super- Arrhcnius 
manner, see Fig. [2] The NEF model behaves as a fragile 
supercooled liquid, as expected Note however that 
the ratio T a /T g w 7 that can be extracted from Fig. [21 
is about three times smaller in fragile liquids such as or- 
thoterphenyl. 

Various fits are also included in Fig. [21 The high tem- 
perature behaviour is well described by a naive mean-field 
approximation plL I23I] . 

t M f ~ c -1 . (6) 

This behaviour breaks down below T Q , where dynam- 
ics is dominated by local fluctuations of the mobility. 
Generalizing to three dimensions the results for the East 
model 0, Il3. l2*4j , we expect that in the limit of small 
temperature, the time scale behaves as 



with A(T) « b/T at low T, where b is a constant. 
Since c ~ e" 1 / 7 , Eq. Q implies an exponential inverse 
temperature squared dependence at low temperatures, 
t ~ e b / T . In fact, an empirical fit using 

A(T) = o + 4 (8) 
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FIG. 2: Points: Mean relaxation time from P(t) = e _1 . 
The lines represent various fits. Dashed line is the high- 
temperature Arrhenius behaviour, r ~ exp (1/T). Full line is 
t(T) ~ exp(A(T)/T) with A(T) = 0.634/T+l.l. The dotted 
line shows that a Vogel-Fulcher form, r ~ exp[2.3/(T — 0.06)], 
fits the low temperature data with an unphysical finite T sin- 
gularity. 
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FIG. 3: Distribution of persistence time •nit), see Eq. Jjjj, for 
temperatures as in Fig.0 decreasing from left to right. 



titatively reproduced by the NEF model. 

B. Distribution of relaxation times 



allows one to describe the whole temperature range stud- 
ied, using b = 0.634 and a = 1.1. The supposedly exact 
low-temperature result, A(T) = (31n2T) _1 , works well 
for very large relaxation times. Indeed we find that a 
plot of T 2 logr vs. 1/T nicely converges to the theoret- 
ical value (3m2) _1 (not shown). At low temperatures, 
the a-relaxation in the NEF model therefore follows a 
Bassler law [25| . 

Fig. |21 also shows that the low T data can be fit us- 
ing the Vogel-Fulcher law, t vf ~ exp [A/(T - T )]. The 
Bassler and Vogel-Fulcher curves are hardly distinguish- 
able for over eight orders of magnitude, as has been ob- 
served before when fitting experimental data [25j . In our 
case, it is evident that the finite temperature singularity 
of the Vogel-Fulcher law is completely unphysical. 

As will be discussed below, the fragile behaviour of 
the NEF model results from the existence of a hierarchy 
of lengthscales governing the dynamical behaviour. It is 
therefore not surprising that the relaxation time in this 
model does not follow the Adam-Gibbs relation which is 
argued to link dynamics and thermodynamics of super- 
cooled liquids |2|||23]. For the NEF model, the entropy 
has essentially an Arrhenius behaviour at low temper- 
atures. The Adam-Gibbs relation thus predicts the re- 
laxation time to grow as a double exponential of the in- 
verse temperature. This vastly overestimates the growth 
of r, which follows the Bassler law, Eqs. I|7I8|I . In the 
same vein, naively extracting a Kauzmann temperature 
by linearly extrapolating the entropy to zero in the same 
temperature range where the Vogel-Fulcher law seems to 
apply (see Fig. EJ) , gives a value Tk ~ 0.16, which is very 
different from the Vogel-Fulcher temperature To as 0.06 
found above. Although Tq and Tk are ill-defined temper- 
atures in our case, experimentally they are often found 
to be close in fragile liquids, a feature which is not quan- 



The mean relaxation time r(T) captures only in part 
the relaxation behaviour of the model. We consider in 
this subsection the distribution of relaxation times, 7r(t), 
related to the mean persistence function via |2l| 



P(t) 



dt'ir{t'). 



(9) 



These distributions are shown in Fig. [3J 

As mentioned above, the long-time tail of the per- 
sistence function is well-described by a stretched expo- 
nential form, which implies the following long-time be- 
haviour for the distribution of relaxation times, 



■ exp 



(10) 



valid for when t > t. As for the East model |4 l28ll29| . 
we find that the stretching exponent decreases with T 
in the regime T < T a , from its high temperature value 
0=1. 

As can be guessed from Figs. ^ and using r(T) and 
0(T) as unique fitting parameters does not allow for a 
satisfactory description of the whole decay of the persis- 
tence function and distribution of relaxation times. This 
becomes evident when data are presented in an alterna- 
tive way. Following the experimental literature, we show 
the data obtained from the distribution of time scales in 
a frequency representation via |30j 



x"M = Im 



7r(log(r)) 



1 



1 



dlogr, (11) 



as is often done with dielectric susceptibility measure- 
ments. This will also allow us to make quantitative com- 
parisons to experiments below. The results for x"(w) are 
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FIG. 4: Imaginary part of the dynamic susceptibility denned 
in Eq. li lt , for temperatures as in Fig. decreasing from 
right to left. For T = 0.18, we also show as a dashed line 
the Fourier transform of the stretched exponential fit to the 
long time tail of the persistence function, revealing 'additional 
processes' on the 'high-frequency flank' of the a-relaxation. 



displayed in Fig. 0] At temperature T = 0.18 we show 
on top of the data the Fourier transform of the stretched 
exponential fit to the long time tail of P(t). The fit is 
clearly incorrect in this representation, although it would 
look acceptable in most of a plot such as in Fig. ^ I n 
fact, deviations appear when cut > 1, corresponding to 
t < t in the time domain. This reveals the presence of 
'additional processes' on the 'high-frequency flank' of the 
a-relaxation. 

This feature is obviously reminiscent of the 'high- 
frequency', or 'Nagel' wing extensively studied by di- 
electric spectroscopy in several materials |3l|. The wing 
is usually observed in fragile glass-formers, but its pre- 
cise nature has not been fully elucidated. Despite ini- 
tial claims of certain universal behaviour of the high- 
frequency wing, more recent investigations seem to favour 
the interpretation that the phenomenon does not obey 
universal scaling (3l| . This is compounded by the obser- 
vation that the amplitude of the phenomenon also seems 
to depend on the technique used to study it [32] ]. 

In our case, the physical interpretation of the wing is 
very clear. The relaxation of the NEF model proceeds 
in a hierarchical manner, just as in the one-dimensional 
case 0, . In order to relax a domain of immobile regions 
of size £, a region of size 1/2 must first be relaxed, which 
itself necessitates the relaxation of a domain of size £/4, 
etc., down to the smallest size £ = 1. This implies an 
energy cost AE(£) ~ \n£, from which the Bassler form, 
Eq.JSJ, and stretched exponential decay, Eq. QlOfl. fol- 
low |6j. Interestingly, the hierarchy also shows up in the 
distributions of Figs. |31 and 0] at time scales shorter than 
the a-relaxation in a manner reminiscent of the experi- 
mental finding of additional short-time processes. Due to 
the underlying lattice structure of the model, the hierar- 
chy is discrete rather than continuous, as can be indeed 
observed in Fig. Notice also that there is no 'fast' 



process at the microscopic time scale, u> ~ 1, in our sys- 
tem. This is a consequence of coarse-graining by which 
molecular vibrations are removed. 

These results show that the pattern of dynamical re- 
laxation of the NEF model is quantitatively accurate on 
a wide frequency range. This will be shown explicitly 
in Sec. IVI where dielectric susceptibility data taken close 
to the experimental glass transition are compared to the 
NEF model results. We recall that in the strong case 
the large-time decay was found to be purely exponential, 
while short-time processes had a much smaller magni- 
tude |l4j . This prediction is however hard to confirm 
experimentally since strong liquids are not easy to study 
by dielectric spectroscopy |33j. 

Physically, our results also suggest that the high- 
frequency wing observed in fragile glass-formers is an in- 
trinsic feature of the a-relaxation linked in a direct way 
to their dynamical hierarchical structure which is also at 
the origin of stretched relaxation and fragile behaviour. 

When the relaxation time is moderate, an even more 
complex behaviour is observed due to the influence of 
clusters of defects. Such dynamic objects are irrelevant 
at low temperatures, but can influence the short-time dy- 
namics just below T , as discussed in detail in Ref. |21| . 
These clusters were in particular shown to be responsible 
for the temperature behaviour of a number of quantities 
discussed in some numerical works. Clusters also pro- 
duce additional short-time processes visible in the distri- 
butions 7r(t) when temperature is not too low |2l|. In- 
terestingly in the present context, these patterns closely 
resemble the ones observed experimentally in mildly su- 
percooled liquids and usually explained in terms of mode- 
coupling theory 34] . This will be demonstrated in Sec. IVI 
where data recently fitted via mode-coupling theory |35| 
are also fitted using the NEF model. 



IV. LENGTH SCALES 
A. Dynamic heterogeneity 

The growth of timescales in the NEF model is accom- 
panied by growing spatial correlations as the system ap- 
proaches its critical point at T = 0. These correlations 
are purely dynamical in origin, and give rise to dynamic 
heterogeneity [HE!- Figures |3] and serve to illustrate 
this phenomenon, as we now explain. 

We quantify the local dynamics via the local persis- 
tence function Pi{t). For a given temperature we run 
the dynamics for a time t* , such that P(t*) = 1/2, 
meaning that half of the sites have flipped at least once. 
We colour white persistent (immobile) spins, for which 
Pi(t*) = 1, and black transient (currently or previously 
mobile) spins, for which Pi(t*) = 0. Figure [5] shows 
the local persistence function for the NEF model at two 
different temperatures, T G = 1.0 and T = 0.15 <C T G . 
Clearly, the low-temperature dynamics is spatially het- 
erogeneous, and the spatial correlations of the local dy- 



FIG. 5: Spatial distribution of the local persistence at time 
t* such that P(t*) — 1/2 (i.e., 50% of sites, shown in black, 
have flipped by time t*) at temperatures T = 1.0 (top) and 
T — 0.15 (bottom) for a system size L = 40 in both cases. 
The appearance of dynamic critical fluctuations when T — > 
is evident. 



namics grow as T is decreased. The 'critical' nature of 
dynamic clusters is apparent: the pictures are reminis- 
cent of the spatial fluctuations of an order parameter 
close to a continuous phase transition, such as the mag- 
netization of an Ising model near criticality. In our case, 
the order parameter is a dynamic object, the persistence 
function, and the critical fluctuations are purely dynam- 
ical in origin [Isj . 

It is interesting to note that these figures are qualita- 
tively different from the ones obtained in the strong case 
where dynamic facilitation is isotropic. One can clearly 
distinguish in Fig. [S] the North, East and Front direc- 
tions of facilitation, implying that wandering of excita- 
tions in the other three directions is forbidden. Domains 
of Fig. [S] appear much less rough than the ones obtained 
in the isotropic case In that sense, increasing the 

fragility is similar to increasing the 'surface tension' of the 
dynamic domains observed in Fig. |SJ The same observa- 
tion applies to an even more fragile system, the two-spin 
facilitated FA model in two dimensions, where the cor- 



FIG. 6: Visualization of 5% of the sites that relax faster (top) 
and slower (bottom) in a given run at T = 0.18 for a system 
size L = 30. Fast/slow sites are painted as spheres to highlight 
the similarity between heterogeneity in lattice models to that 
in atomistic models. The fast clusters in the top figure are 
less compact than the slow ones in the bottom figure as a 
consequence of relaxation via point like excitations. See the 
text for details. 



responding domains resemble a polydisperse assembly of 
squares pj- 

The qualitative observations of dynamic heterogeneity 
performed in numerical or experimental works can also 
be made in the present coarse-grained model. We show 
in Fig.HjJthe analog of spatial clustering of subsets of fast 
and slow sites. To build these snapshots, we represent 
a given percentage, 5%, of the sites that relax faster or 
slower, i.e. we show those sites for which the persistence 
time is among the 5% smaller or larger in a randomly 
chosen run at temperature T — 0.18. 

One observes that the fastest sites are not randomly lo- 
cated in space, but clustered in 'non-compact' or 'stringy' 
objects, similar to those observed in simulations and ex- 
periments H^, IU |38|. The shape of these objects is a 
consequence of the existence of point defects of mobil- 
ity. When a defect moves, it induces those sites along its 
trajectory to relax, so leaving in its wake a string of fast 
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sites. 

In Fig. we also show 5% of the sites which are slow- 
est, using the same simulation temperature as before. A 
more compact structure is seen. This is again the con- 
sequence of the relaxation via point defects of mobility. 
The slowest sites belong to regions of space devoid of 
defects which take then a very long time to be visited 
by defects. These large domains are thus slowly relaxed. 
It is the bulk of these slow domains that is observed in 
Fig. Note finally that at large times, the distribution 
of slow cluster sizes seems very wide, since some isolated 
sites which have been not been visited by defects coexist 
with the very large domains discussed above. 

Since our comments on Fig. E| are mainly qualitative, 
it should come as no surprise that snapshots built in this 
fashion in the isotropic case are very similar |14| . 

Finally, it is interesting to compare these figures to 
those in previous publications 0, l2l| . in which space- 
time diagrams of one dimensional kinetically constrained 
models were presented. There, spatio-temporal 'bubbles' 
of immobile regions, bounded by diffusing point defects 
were presented. In the present model these bubbles be- 
come (3 + l)-dimensional objects, and Fig.^their three- 
dimensional spatial projections of trajectories of various 
time extensions. 



B. Spatial correlations 

We now quantify these qualitative observations via ap- 
propriate statistical correlators. We measure spatial cor- 
relations of the local dynamics via the following spatial 
correlator of the local persistence function, 



CM,T) 



1 



Nf(t) 



N 

^[mn+rW-p'it)], (12) 



where the function f(t) = P(t) — P 2 {t) in the denomi- 
nator ensures the normalization C(r = 0,t,T) — 1. Al- 
ternatively, one can take the Fourier transform of 1)12(1. 
giving the corresponding structure factor of the dynamic 
heterogeneity, 



1 N 

S(q,t,T) = Jjjjjsj ^2 [{Pk(t)Pi(t))-P 2 (t) 



k,l=l 



Jq(k-l) 



for wavevectors defined in the first Brillouin zone of the 
cubic lattice. 

Finally, the zero wavevector limit of S(q, t, T) defines 
a dynamic 'four-point' susceptibility, x(t,T) = S(q = 
0,t,T), which can be rewritten as the normalized vari- 
ance of the (unaveraged) persistence function, p{t) = 

X (t,T) = [(/(£)> -M*)> 2 ]- (13) 

It should be noted that normalizations also ensure that 
x{t, T) remains finite in the thermodynamic limit, except 




FIG. 7: Dynamic 'four-point' susceptibility, Eq. dl:>> . at tem- 
peratures T = 1.0, 0.6, 0.4, 0.3, 0.23, 0.2, 0.18, and 0.16 (from 
left to right). 



at a dynamic critical point such as the ones discussed in 

Refs. mm. 

Figure shows the time dependence of the suscepti- 
bility l|13fl for various temperatures. The behaviour of 
X is similar to that observed in atomistic simulations of 
supercooled liquids The susceptibility develops at 
low temperature a peak whose amplitude increases, and 
whose position shifts to larger times as T decreases. As 
expected, the location of the peak scales with the relax- 
ation time t(T), indicating that dynamical trajectories 
are maximally heterogeneous when t « t(T). 

In Figure [S] we show the correlator C(r,t,T) and the 
structure factor S(q, t, T) for different temperatures and 
fixed times t = t(T) where dynamic heterogeneity is 
maximal. These correlation functions confirm, as sug- 
gested by Fig. \5\ that a dynamic length scale associated 
with spatial correlations of mobility develops and grows 
as T decreases: spatial correlations decay more slowly 
with distance as T decreases, while a peak in the struc- 
ture develops and grows at q = 0. 

It is possible to extract numerically the value of the 
corresponding dynamic length scale, £(T), at each tem- 
perature. To do so, we study in detail the shape of the 
correlation functions shown in Fig. [S] As for standard 
critical phenomena, we find that the dynamic structure 
factor can be rescaled according to 



%,T,T)~ X (r,T)5(gO 



(14) 



where the scaling function S(x) has the following asymp- 
totic behaviours in the scaling regime, £ 3> 1: 



S{x — > 0) ~ const, 
S{x —* oo) 



■2- v 



(15) 
(16) 



Both the susceptibility \ and the dynamic length scale 
£ estimated at time t = r behave as power laws of the 
defect concentration, 



(17) 
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FIG. 8: Top: Spatial correlator of dynamic heterogeneity at 
the relaxation time t(T). Temperatures as in Fig. Q decrease 
from left to right. Bottom: corresponding structure factor. 
Temperatures as in Fig. decrease from bottom to top. 




FIG. 9: Top: Spatial correlator rescaled with the form C(r) ~ 
C(r/£), £ ~ c -1 / 3 . Bottom: Rescaled structure factor with 
X = S(q = 0) ~ c" 1 and £ ~ c~ ly/3 . Full line interpolates 
between S(g — » 0) ~ const and S^g — » oo) ~ fc -3 ' 58 . 



These relations imply that the exponents 7 and v should 
be numerically accessible by adjusting their values so that 
a plot of c 1 S versus qc~ v is independent of temperature. 
We show such a plot in Fig.03 and we find that the values 
7 « 1 and v « 1/3 lead to a good collapse of the data. 
The exponent 7 can be independently and more directly 
estimated from Fig. [7| by measuring the height of the 
maximum of the susceptibility for various concentrations. 
The result is also well fitted to the power law c 7 with 
7 « 1. We find also that the scaling function S(x) is 
well-described by an empirical form S(x) = l/(l + x 2 ~ v ), 
consistent with Eq. (|16|) . Thus we can determine the 
value of the 'anomalous' exponent, 77. We find that the 
theoretically expected rj = — In 3/ In 2 « —1.58 describes 
the data very well [6j. The fact that this exponent is 
much more negative than in the strong case of the FA 
model quantifies our observation that domains of Fig. 
were much less rough in the NEF model than in the FA 
model, i.e., that dynamic domains of fragile systems have 
smoother boundaries than those of strong ones. 

The spatial correlator C(r, r, T) also obeys scale invari- 
ance, and we find that 

C(r,T,T)~c(j,r) , (18) 



with the asymptotic behaviour C(x — > 0) ~ a: -0 ' 58 , as 
a result of a generalized Porod's law. The scaling be- 
haviour l|18[) is presented in Fig. which also confirms 
the temperature dependence of the correlation length, 
£ ~ c -1 / 3 12] . It is interesting to note that the relation 
between susceptibility and length is very different in the 
strong and fragile cases, since we find that x ~ £ 3 m the 
NEF model, while the scaling is closer to x ~ £ 2 m the 
FA case 0. 



V. COMPARISON TO EXPERIMENTAL DATA 

In this section, we use the numerical results obtained in 
previous sections to fit experimental data with the NEF 
model. In doing so there are several points that need to 
be considered. 

First, the model we use is meant to be a description of 
supercooled liquids which is coarse-grained both in time 
and space, and lives on a lattice. We are thus dealing 
with discrete, rather than continuous, spatial degrees of 
freedom, and the very short-time dynamics of the liquid 
is removed. By construction, this produces discrepancies 
between real data on liquids and numerical NEF model 
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FIG. 10: Optical Kerr Effect experimental for supercooled 
liquid 2-biphenylmethanol (full lines from Ref. |4l|) at tem- 
peratures T = 359 K, 327 K, 319 K, 311 K, 303 K, 291 K 
(from left to right). Dashed lines are corresponding data for 
the NEF model at temperatures T = 0.4, 0.34, 0.32, 0.3, 0.28 
and 0.25 (from left to right). Time is counted in picoseconds 
in the experiments, in Monte Carlo steps in simulations. 



data, especially for short times and small lengths. This 
is a small price to pay in such an approach given the 
large number of features that can still be satisfactorily 
accounted for with the NEF model. 

Second, we have some freedom on how to relate real 
experimental timescales and Monte Carlo steps in the 
simulation. This will be done empirically, and we find 
that the expected equivalence, 1 MC step « few ps [l^ j. 
works well, independently of the temperature. Explicitly, 
we found that 1 MC <E [1, 10] ps for the whole range of 
experimental data we have considered 0, . 

Third, one must adjust the temperature given in 
Kelvin in experiments to the adimensional T of the sim- 
ulations. This amounts to fitting the value of an energy 
scale, J, which should appear on dimensional grounds 
in front of the Hamiltonian In principle, J could 

be fixed independently by fitting, say, viscosity data of 
a given liquid before using the corresponding tempera- 
tures to fit more detailed dynamic data in what would 
be a zero-parameter fitting procedure. We have done the 
correspondence in a less constrained manner, adjusting 
T in the simulation to give a good fit to the data. Very 
satisfactorily, though, we end up with a correspondence 
between numerical and experimental temperatures which 
is well described by linear relations, as it should be at low 
temperatures .12]- 

In Fig. we show the result of this fitting procedure 
as applied to recent data measured using the so-called 
optical Kerr effect EH- This technique has several 
advantages. It extends for over five orders of magnitude 
in time. The quantity measured is the derivative of a 
time correlation function, and therefore the analog of the 
distribution of time scales, 7r(i), discussed above. And vi- 
brations, which are neglected in our approach, affect very 
little the measured decay in the timescales of interest. 



In Ref. |4l|, the experimental results were fitted using 
the empirical form 

7r(t) = [pt- 1+c + dt^ 1 } exp(-Vr), (19) 

with the values i«0.8 - 0.85 and C » 0. - 0.2 con- 
sistently found in different liquids. The l/t behaviour 
of 7r(t) obtained at short times when C = implies a 
logi behaviour of the time correlator. This was taken as 
a challenging result to mode-coupling theory which does 
not naturally predict such patterns, as can be seen when 
asymptotic analytic results are considered [34|]. These 
data were however recently virtually perfectly fitted us- 
ing a standard schematic version of the mode-coupling 
theory, therefore nullifying the criticisms raised in the 
experimental work j3j| . Note that several fitting param- 
eters are allowed by the mode-coupling fits, which are 
performed using a two-correlator schematic model: cou- 
pling constant, distance to the dynamic singularity and 
various numerical factors adapting theoretical time scales 
to the experimental ones. By contrast, we make use of 
just one free parameter in Fig. 

The 'nearly logarithmic decay' of correlations has a 
simple explanation in our approach. At modest temper- 
atures, T < T Q , there is a coexistence of isolated ex- 
citations, responsible for 'slow processes' and rapidly 
relaxing clusters of excitations; see Ref. |2J] for a de- 
tailed discussion of the related temperature crossovers. 
This coexistence manifests as a nearly flat distribution 
7r(logt) over a wide range of timescales. This translates 
into a l/t behaviour of n(t), and logarithmic decay of 
P(t), as observed in Fig. This explanation, more- 
over, predicts that the effective exponent C appearing in 
Eq. I|19|) should acquire a slight temperature dependence, 
and change from C < 0, when fast processes dominate 
closer to T a , to C > at lower temperature when slow 
processes become dominant. This is precisely what is 
found in experiments • This subtlety is also accounted 
for by mode-coupling theory |35l ]. 

The only real discrepancy between fits and data is vis- 
ible at very short time, as was anticipated, but the over- 
all quantitative agreement is very good. As already sug- 
gested by the qualitative analysis of Ref. [2lJ , our theoret- 
ical approach can be used even far above the experimen- 
tal glass transition. Mode-coupling theory is thought to 
be applicable in this regime, under the assumption that 
the dynamics at modest supercooling is different from 
that near T g . As in Ref. |2l|, our results here suggest 
that this may not be the case. 

One of the advantages of the present approach on 
mode-coupling theory is that it does not produce an un- 
physical singularity at a temperature above T g , and it 
can therefore be used down to very low temperatures and 
large relaxation times. In Fig. ^2 we use the NEF model 
to fit dielectric data taken on Salol [3l| . As before, the fits 
are done using a single free parameter. Fig.llllshows that 
the overall agreement is again very good. Note in partic- 
ular that the high-frequency wing is correctly accounted 
for by the (discrete) hierarchy of dynamic lengthscales 
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FIG. 11: Dielectric susceptibility data for Salol obtained from 
the authors of Refs. JJi at temperatures T = 2557f, 243K, 
233if, 225JS", and 219 K (from right to left). Full lines are 
corresponding data for the NEF model at temperatures T = 
0.385, 0.26, 0.198, 0.168, and 0.15 (from left to right). Time 
is counted in units of 10 ps in the experiments, and in Monte 
Carlo units in the simulations. 



discussed in the previous sections. 

Again, discrepancies due to coarse-graining are evi- 
dent: absence of short-time processes, and discreteness 
of the hierarchy of time scales. Discrete scales are par- 
ticularly evident in the Nagel wing, where numerical 
data only produce the 'skeleton' of the wing instead of 
a smooth curve. Note that the same experimental data 
were fitted in Ref. 42] using a frustration limited do- 
main scaling picture of the glass transition: seven fit- 
ting parameters were used there to obtain a satisfactory 
fit. Such data are usually fitted in experimental papers 
by empirical forms involving again several free parame- 
ters [Ig. 

Consider, finally, the scaling of lengths and times. Fig- 
ure 1121 shows dynamic scaling in three different model 
systems: the three-dimensional FA model (indicated as 
'strong'; data from Refs. 0]), the present NEF model 
data (indicated as 'fragile'), and the Kob- Andersen 
Lennard-Jones binary mixture |43T ] ('LJ' in the figure; 
data from Ref. The figure shows that in the fragile 

case the growth of dynamic lengthscales is much less pro- 
nounced than in the strong one. This is one of the central 
predictions of the dynamic facilitation approach [l^. It 
is a consequence of the temperature dependent dynamic 
exponent of East-like models such as the NEF model, in 
contrast to strong systems where dynamic lengths go as 
a fixed power of the relaxation time. Note that a figure 
qualitatively similar to Fig.^|would be obtained for the 
four-point susceptibility as a function of relaxation time. 

This slow growth of lengthscales in the NEF model is 
also compatible with the modest size of heterogeneities 
as suggested by experiments near T g [s|. Interestingly, 
the LJ mixture appears in this representation as a fairly 
strong system, despite the common assumption that it is 
a good model for fragile liquids. The theoretical and nu- 




FIG. 12: Dynamic scaling of timescales versus lengthscales 
in the strong three-dimensional FA model |l3j|. the fragile 
three-dimensional NEF model, and three-dimensional binary 
Lennard-Jones mixture [l3T| . 



merical findings reported in Fig. 1 121 offer a solution to the 
paradox that experimentally measured length scales are 
very small when compared to extrapolations performed 
from numerical works |8l III . 



VI. CONCLUSIONS 

In summary, we have presented an extensive numeri- 
cal study of the North-or-East-or- Front model, the three- 
dimensional generalization of the East facilitated spin 
model. The NEF model, whose dynamical rules have an 
externally imposed asymmetry, is in the same universal- 
ity class of the more physical arrow model of Ref. ^| m 
the limit of maximal directional persistence, and there- 
fore models the behaviour of fragile, super- Arrhenius, 
glass formers. 

We have characterized the equilibrium dynamics of the 
NEF model through measurements of the cv-relaxation 
timescales, distributions of relaxation times, and the 
growth of dynamic correlation lengths and four-point 
susceptibilities. We find that r Q follows a Bassler law 
of exponential-inverse-squared temperature dependence, 
and that the dynamic correlation length grows as £ ~ c~ u 
with v ks 1/3, where c is the equilibrium concentration 
of excitations, while the dynamic susceptibility grows as 
X ~ c~ 7 with 7 « 1. These results confirm some of the 
expectations of Ref. [12| for the arrow model, in partic- 
ular the quasi one-dimensional nature of the dynamics 
in all dimensions. The decay of the four-point structure 
factor, S(k) « A: -2-158 , is a direct generalization of the 
result for the ID East model @, and an indication of 
hierarchical dynamics in the NEF model. These two fea- 
tures, persistence of directionality and hierarchical dy- 
namics, give rise to fragile behaviour |l2j |. 

We have shown how the NEF model can be used to ra- 
tionalize experimental observations over a wide range of 
temperatures and timescales. We have fitted correlation 
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data in the time domain, obtained in the mildly super- 
cooled regime by optical Kerr effect ^l|> with the NEF 
model predictions. This provides an interpretation of the 
quasi logarithmic time dependence of correlations based 
on the coexistence of isolated excitations and clusters 
of excitations, which occurs when the system is crossing 
over from a regime of homogeneous dynamics at high T 
to one of heterogeneous dynamics at low T |2l| . We have 
shown how the NEF model can also be used to describe 
dielectric susceptibility data, in the frequency represen- 
tation, measured close to the glass transition by dielectric 
spectroscopy. Interestingly, the NEF model successfully 
accounts for the excess high-frequency or Nagel wing [31| • 
The resul ts p resented in this paper add to the list 
[3 EU l2fi l44| of phenomenological observations that 
can be quantitatively understood within the dynamic fa- 



cilitation approach. 
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